PROGRAM main
IMPLICIT NONE

!!!!!!!!!!!!!!!!!!!!
! DEFINE VARIABLES !
!!!!!!!!!!!!!!!!!!!!

INTEGER :: i1, i2, i3, i4, i5, i6, i7, i8, i9, i5F, iy, i_inc
!it=1: before shock, it=2: after shock
!group 1: low robot risk/ group 2:medium robot risk / group 3: high robot risk 
INTEGER :: aux2, ig, it 
INTEGER :: tb, t1, tr, td, tn, tw, tn1, t, t2, per, nper !nper: number of periods after robotization shock ("the robotization phase")
integer :: ng, nt !ng: number of groups, nt: "periods", before and after robotization
parameter (tb=22, tr=66, td=100, tw=52-tb+1, tn=td-tb+1, nper=8, ng=3, nt=2) 
INTEGER :: na, nlw, n, nc, ninc
!parameter (na=15, nlw=51, ninc=9, n=3, nc=16)
parameter (na=21, nlw=51, ninc=15, n=3, nc=16)
INTEGER, DIMENSION(1) :: pt
REAL :: maxlw=50.0, minlw=0.0001 
REAL :: l_maxlw, l_minlw, steplw
real :: max_inc = 5.0, min_inc=0.5, stepinc
REAL :: stepc, maxc, minc, mpc
real, dimension(ng,nt) :: ret_fac0
real, dimension(ng,nt) :: aa, b1, b2, b3, smay, smav   
real, dimension(ng,nt) ::  p_unemp
REAL :: corr_v=0.0, corr_y=0.0 
real :: gs_L = 4375.0/30375.0, gs_M= 16375.0/30375.0, gs_H= 9625.0/30375.0 !group sizes
REAL :: rho_L=6.0, delta_L=0.994, psi_L=0.4 
REAL :: rho_M=6.0, delta_M=0.975, psi_M=0.275  
REAL :: rho_H=6.0, delta_H=0.98, psi_H=0.2
real :: rho, delta, psi, theta, psi_1, psi_2
real :: ct, caux 
real :: py_unemp=-0.15 !permanent income shock following an unemployment state.

real :: Radj=0.01, auxR !adjustment to realized returns in the post-sims (to match the realized returns in the data)
real :: W0=0.0 !initial endowment of wealth, as a fraction of initial income 
REAL :: r=1.01, mu=0.051 !giving an unconditional equity premium of 4% (after disaster risk)
real :: p_lowr=0.02, lowr=0.515 !from Fagereng et al JF
real :: sigr=0.3
!REAL :: F=0.04, Fy=0.0 !group 1
REAL :: F=0.005, Fy=0.0 !group 2
real :: yu1=0.8, yu1_cap=0.1496, yu2_cap=0.064 !unemp. income is divided by 100K, just like labor income below
real :: pct_ins=0.7, unemp_dur=0.358333 !% of workers with unemployment insurance and average unemployment duration
REAL :: net_cash, partcount, income
double precision :: auxVV, u, int_V, Vp
INTEGER :: i_net_cash, i_net_cash2, ic1, ic2, auxpu
REAL :: ttc, sav, rand_u, rand_u2, auxgy, survprob
REAL, DIMENSION(tn-1,1) :: survprob0=0.0, survprob1=0.0
REAL, DIMENSION(n,1) :: grid, weig, ones_n_1=1.0, grid2
REAL, DIMENSION(n+1,1) :: gret, nweig2 
REAL, DIMENSION(n,n+1) :: yp !plus 1 in the second dimension is for the return disaster state
REAL, DIMENSION(n,1) :: yh, w1  !unemployment state is captured separately (hence "just" n) 
REAL, DIMENSION(n,n+1,n+1) :: nweig1
REAL, DIMENSION(n+1,n+1) :: nweig1b
REAL, DIMENSION(ninc,n,n+1,tr-tb) :: inc_t1, tty_t1
INTEGER, DIMENSION(ninc,n,n+1,tr-tb) :: inc_it1
REAL, DIMENSION(ninc,tr-tb) :: inc_t1_u, tty_t1_u
INTEGER, DIMENSION(ninc,tr-tb) :: inc_it1_u
REAL, DIMENSION(tr-tb+1,ng,nt) :: f_y=0.0
REAL, DIMENSION(nlw,1) :: glw, lglw
REAL, DIMENSION(na,1) :: ga
REAL, DIMENSION(ninc,1) :: ginc
REAL, DIMENSION(na,n+1) :: riskret
REAL, DIMENSION(nc,1) :: gc
real, dimension(ng,1) :: gs !group sizes (%)
double precision, DIMENSION(na,nc) :: auxV
double precision, DIMENSION(nc,1) :: auxV2
double precision, DIMENSION(na*nc,1) :: vec_V
real, DIMENSION(nlw,ninc,n+1,ng,nt) :: secd   
REAL, DIMENSION(nlw,tn,ninc,n+1,ng,nt) :: C=0.0, V, A=1.0 

INTEGER :: w, nsim, auxP, indYT1, indYT2
REAL :: cash, aux1, count, simTY1, simTY2, simPY_1, simPY_2
REAL, DIMENSION(1,1) :: eps_r, eps_y
!parameter (nsim=5000)
parameter (nsim=10000)
!parameter (nsim=1)
REAL, DIMENSION(nsim,1) :: ones_nsim_1=1.0
REAL, DIMENSION(tw,1) :: meanYs, dweig=0.0
REAL, DIMENSION(tw+nper,ng) :: meanY, meanY_R
REAL, DIMENSION(tw,ng) :: meanC, meanW, meanS, meanB, meanWY, meanalpha, meanP, meanBpart
REAL, DIMENSION(tw,nsim) :: simC, simW=0.0, simA, simS, simB, simW_Y=0.0, simP=0.0, simBpart, simW_R=0.0 
REAL, DIMENSION(tw,ng) :: meanW_R, meanW_R2

REAL, DIMENSION(tw,1) :: meanYaux 
REAL, DIMENSION(6,1) :: sim_out

real, dimension(ng+1,1) :: mWY=0.0, mP=0.0, malpha=0.0 !+1 is for global average
real, dimension(ng+1,1) :: mW=0.0, mW_R=0.0, mW_R2=0.0 !+1 is for global average
REAL, DIMENSION(tw+nper,nsim,ng) :: simPY, simPY2, simY, tty
INTEGER, DIMENSION(tw+nper,nsim,ng) :: indy, indYT
REAL, DIMENSION(tn,nsim) :: simR 
REAL, DIMENSION(nper,tw,nsim,ng) :: simY_R, simPY_R, simPY2_R, tty_R
REAL, DIMENSION(tw,nsim,ng) :: auxmyr
INTEGER, DIMENSION(nper,tw,nsim,ng) :: indy_R, indYT_R

real, dimension(nper,nsim) :: simW_aux=0.0
REAL :: simP_aux=0.0, simBpart_aux=0.0, simC_aux=0.0, simA_aux=0.0, simS_aux=0.0, simB_aux=0.0, simW_Y_aux



!!!!!!!!!!!!!!!
! GROUP SIZES !
!!!!!!!!!!!!!!!
gs(1,1) = gs_L
gs(2,1) = gs_M
gs(3,1) = gs_H

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! PARAMETERS OF LABOR INCOME PROCESS !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!pre-robotization period
!low robot risk
aa(1,1)= -225515.63950
b1(1,1)=21716.43581
b2(1,1)=-3219.49650/10
b3(1,1)= 193.65901/100
smay(1,1)=0.0864
smav(1,1)=0.0801
p_unemp(1,1)=0.088 
ret_fac0(1,1) = 0.5586 

!medium
aa(2,1)= -208473.29755
b1(2,1)= 19841.56108
b2(2,1)= -2514.56761/10
b3(2,1)=  150.46336/100
smay(2,1)=0.0871
smav(2,1)=0.0701
p_unemp(2,1)=0.0829
ret_fac0(2,1) = 0.5025 

!high robot risk
aa(3,1)= -129141.48317
b1(3,1)= 13862.27457
b2(3,1)= -1259.19925/10
b3(3,1)=  67.73289/100
smay(3,1)=0.0757
smav(3,1)=0.0615
p_unemp(3,1)=0.068
ret_fac0(3,1) = 0.4691 

!post-robotization period
!low robot risk
aa(1,2)= -290334.41946
b1(1,2)=24498.60538
b2(1,2)=-3145.08207/10
b3(1,2)=159.89917/100 
smay(1,2)=0.1130
smav(1,2)=0.0637
p_unemp(1,2)=0.0702 
ret_fac0(1,2) = 0.5873 

!medium robot risk
aa(2,2)= -256002.39436
b1(2,2)=  23692.50904
b2(2,2)= -2854.63861/10
b3(2,2)= 135.59702/100
smay(2,2)=0.1094
smav(2,2)=0.0664
p_unemp(2,2)=0.0645
ret_fac0(2,2) = 0.5638 

!high robot risk
aa(3,2)= -272652.39247
b1(3,2)=  26189.74566
b2(3,2)= -3641.31613/10
b3(3,2)= 187.17217/100
smay(3,2)=0.0963
smav(3,2)=0.0696
p_unemp(3,2)=0.0665
ret_fac0(3,2) = 0.5725 


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! APPROXIMATION TO NORMAL DISTRIBUTION  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
grid(1,1) = -1.73205080756887
grid(2,1) = 0.0
grid(3,1) = 1.73205080756887

weig(1,1) = 0.16666666666667
weig(2,1) = 0.66666666666667
weig(3,1) = 0.16666666666667

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! CONDITIONAL SURVIVAL PROBABILITIES  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
survprob0(1,1)	=	0.9993
survprob0(2,1)	=	0.9993
survprob0(3,1)	=	0.99926
survprob0(4,1)	=	0.99925
survprob0(5,1)	=	0.99934
survprob0(6,1)	=	0.9994
survprob0(7,1)	=	0.99933
survprob0(8,1)	=	0.99928
survprob0(9,1)	=	0.99169
survprob0(10,1)	=	0.98996
survprob0(11,1)	=	0.99147
survprob0(12,1)	=	0.98868
survprob0(13,1)	=	0.98882
survprob0(14,1)	=	0.98809
survprob0(15,1)	=	0.98726
survprob0(16,1)	=	0.98713
survprob0(17,1)	=	0.98379
survprob0(18,1)	=	0.98334
survprob0(19,1)	=	0.98416
survprob0(20,1)	=	0.98109
survprob0(21,1)	=	0.98116
survprob0(22,1)	=	0.9767
survprob0(23,1)	=	0.97687
survprob0(24,1)	=	0.97356
survprob0(25,1)	=	0.96629
survprob0(26,1)	=	0.96681
survprob0(27,1)	=	0.96801
survprob0(28,1)	=	0.96242
survprob0(29,1)	=	0.95647
survprob0(30,1)	=	0.95677
survprob0(31,1)	=	0.95316
survprob0(32,1)	=	0.94835
survprob0(33,1)	=	0.94567
survprob0(34,1)	=	0.93925
survprob0(35,1)	=	0.93805
survprob0(36,1)	=	0.93585
survprob0(37,1)	=	0.92819
survprob0(38,1)	=	0.91962
survprob0(39,1)	=	0.92084
survprob0(40,1)	=	0.91709
survprob0(41,1)	=	0.92032
survprob0(42,1)	=	0.91853
survprob0(43,1)	=	0.91894
survprob0(44,1)	=	0.85533
survprob0(45,1)	=	0.91077
survprob0(46,1)	=	0.90003
survprob0(47,1)	=	0.90699
survprob0(48,1)	=	0.88391
survprob0(49,1)	=	0.78329
survprob0(50,1)	=	0.90577
survprob0(51,1)	=	0.88631
survprob0(52,1)	=	0.89715
survprob0(53,1)	=	0.88443
survprob0(54,1)	=	0.8882
survprob0(55,1)	=	0.86896
survprob0(56,1)	=	0.8806
survprob0(57,1)	=	0.87224
survprob0(58,1)	=	0.85225
survprob0(59,1)	=	0.85839
survprob0(60,1)	=	0.86976
survprob0(61,1)	=	0.84882
survprob0(62,1)	=	0.84447
survprob0(63,1)	=	0.84352
survprob0(64,1)	=	0.79812
survprob0(65,1)	=	0.78191
survprob0(66,1)	=	0.84477
survprob0(67,1)	=	0.75803
survprob0(68,1)	=	0.72327
survprob0(69,1)	=	0.74053
survprob0(70,1)	=	0.68536
survprob0(71,1)	=	0.63742
survprob0(72,1)	=	0.56165
survprob0(73,1)	=	0.57903
survprob0(74,1)	=	0.41948
survprob0(75,1)	=	0.49974
survprob0(76,1)	=	0.45874
survprob0(77,1)	=	0.16423
survprob0(78,1)	=	0.14104


survprob1(1,1)	=	0
survprob1(2,1)	=	0
survprob1(3,1)	=	0
survprob1(4,1)	=	0
survprob1(5,1)	=	0
survprob1(6,1)	=	0
survprob1(7,1)	=	0
survprob1(8,1)	=	0
survprob1(9,1)	=	0.00063
survprob1(10,1)	=	0.00076
survprob1(11,1)	=	0.00064
survprob1(12,1)	=	0.00086
survprob1(13,1)	=	0.00085
survprob1(14,1)	=	0.0009
survprob1(15,1)	=	0.00097
survprob1(16,1)	=	0.00097
survprob1(17,1)	=	0.00123
survprob1(18,1)	=	0.00126
survprob1(19,1)	=	0.00119
survprob1(20,1)	=	0.00143
survprob1(21,1)	=	0.00142
survprob1(22,1)	=	0.00177
survprob1(23,1)	=	0.00174
survprob1(24,1)	=	0.00199
survprob1(25,1)	=	0.00255
survprob1(26,1)	=	0.0025
survprob1(27,1)	=	0.0024
survprob1(28,1)	=	0.00282
survprob1(29,1)	=	0.00328
survprob1(30,1)	=	0.00325
survprob1(31,1)	=	0.00352
survprob1(32,1)	=	0.00386
survprob1(33,1)	=	0.00407
survprob1(34,1)	=	0.00455
survprob1(35,1)	=	0.00463
survprob1(36,1)	=	0.00478
survprob1(37,1)	=	0.00537
survprob1(38,1)	=	0.006
survprob1(39,1)	=	0.00587
survprob1(40,1)	=	0.00614
survprob1(41,1)	=	0.00585
survprob1(42,1)	=	0.00595
survprob1(43,1)	=	0.00586
survprob1(44,1)	=	0.01096
survprob1(45,1)	=	0.00639
survprob1(46,1)	=	0.00724
survprob1(47,1)	=	0.00659
survprob1(48,1)	=	0.00837
survprob1(49,1)	=	0.01667
survprob1(50,1)	=	0.00647
survprob1(51,1)	=	0.00795
survprob1(52,1)	=	0.00683
survprob1(53,1)	=	0.00769
survprob1(54,1)	=	0.00716
survprob1(55,1)	=	0.00852
survprob1(56,1)	=	0.00728
survprob1(57,1)	=	0.00761
survprob1(58,1)	=	0.00892
survprob1(59,1)	=	0.00792
survprob1(60,1)	=	0.00648
survprob1(61,1)	=	0.00762
survprob1(62,1)	=	0.00726
survprob1(63,1)	=	0.00659
survprob1(64,1)	=	0.00953
survprob1(65,1)	=	0.00989
survprob1(66,1)	=	0.00339
survprob1(67,1)	=	0.00945
survprob1(68,1)	=	0.01093
survprob1(69,1)	=	0.00801
survprob1(70,1)	=	0.01085
survprob1(71,1)	=	0.01341
survprob1(72,1)	=	0.01784
survprob1(73,1)	=	0.0145
survprob1(74,1)	=	0.02627
survprob1(75,1)	=	0.01719
survprob1(76,1)	=	0.01815
survprob1(77,1)	=	0.04234
survprob1(78,1)	=	0.04114

!!!!!!!!!!!!!!!!!!!!!!!
! DEMOGRAPHIC WEIGHTS !
!!!!!!!!!!!!!!!!!!!!!!!
OPEN(UNIT=25,FILE='dem_weig.txt',STATUS='old',ACTION='read')
do i1=1,tw
   read(25,*) dweig(i1,1)
end do !ibg
CLOSE(UNIT=25)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! GRIDS FOR THE STATE VARIABLES AND FOR PORTFOLIO RULE !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do i1=1,n
   gret(i1,1) = r+mu+grid(i1,1)*sigr
end do 
gret(n+1,1) = lowr

do i1=1,n
   nweig2(i1,1) = weig(i1,1)*(1-p_lowr)
end do 
nweig2(n+1,1) = p_lowr

DO i1=1,na
   ga(i1,1)=(na-i1)/(na-1.0)
END DO

do i5=1,na
   do i8=1,n+1
      riskret(i5,i8)=r*(1-ga(i5,1))+gret(i8,1)*ga(i5,1)
   end do   !i8
end do    !i5


l_maxlw = log(maxlw)
l_minlw = log(minlw)
steplw = (l_maxlw-l_minlw)/(nlw-1)
DO i1=1,nlw
   lglw(i1,1)=l_minlw+(i1-1.0)*steplw
END DO
DO i1=1,nlw
   glw(i1,1)=exp(lglw(i1,1))
END DO


stepinc = (max_inc-min_inc)/(ninc-1)
DO i1=1,ninc
   ginc(i1,1)=min_inc+(i1-1.0)*stepinc
END DO
write(*,*) ginc


!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATED RETURNS SHOCKS !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,tn
do i1=1,nsim/2
   call randn(1,eps_r(1,1))
   simR(t,i1) = mu + r + sigr*eps_r(1,1)
   simR(t,nsim/2+i1) = mu + r - sigr*eps_r(1,1)
end do   
do i1=1,nsim
   call random_number(rand_u)
   if (rand_u<p_lowr) then
      simR(t,i1) = lowr
   end if    
end do
end do   

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DETERMINISTIC INCOME PROFILE !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do it=1,nt
do ig=1,ng


DO i1=tb,tr
   f_y(i1-tb+1,ig,it) = aa(ig,it)+b1(ig,it)*i1+b2(ig,it)*i1**2+b3(ig,it)*i1**3
   f_y(i1-tb+1,ig,it) = f_y(i1-tb+1,ig,it)/100000   
END DO

end do
end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATED INCOME SHOCKS !
!!!!!!!!!!!!!!!!!!!!!!!!!!!
!pre-robotization
it=1

do ig=1,ng
!cumulative densitiy function for the income shocks (with the top mass "empty" corresponding to the unemployment state)
w1(1,1) = weig(1,1)*(1-p_unemp(ig,it))
do i7=2,n
   w1(i7,1) = w1(i7-1,1)+weig(i7,1)*(1-p_unemp(ig,it))
end do   
end do

do i1=1,nsim/2
do i2=1,tw+nper !need to simulate additional nper for the counterfactual exercise
call random_number(rand_u) !before ig loop: same shocks for all groups
call random_number(rand_u2) !before ig loop: same shocks for all groups
do ig=1,ng
   if (i2.eq.1) then
      call simYP_f(rand_u2,weig(:,1),grid(:,1),n,smav(ig,it),simPY_1,simPY_2)
      simPY(1,i1,ig) = simPY_1
      simPY(1,nsim/2+i1,ig) = simPY_2
   else
      if (indYT(i2-1,i1,ig).eq.4) then
         simPY_1 = py_unemp
         simPY_2 = py_unemp
      else
         call simYP_f(rand_u2,weig(:,1),grid(:,1),n,smav(ig,it),simPY_1,simPY_2)
      end if
      simPY(i2,i1,ig) = simPY_1+simPY(i2-1,i1,ig)
      simPY(i2,nsim/2+i1,ig) = simPY_2+simPY(i2-1,nsim/2+i1,ig)
   end if 
   simPY2(i2,i1,ig) = exp(simPY(i2,i1,ig))*f_y(i2,ig,it)
   simPY2(i2,nsim/2+i1,ig) = exp(simPY(i2,nsim/2+i1,ig))*f_y(i2,ig,it)
    !for the transitory shocks we consider the unemployment and non-unemployment states separately
   if (rand_u<w1(3,1)) then
      call simY_f(rand_u,w1(:,1),grid(:,1),n+1,smay(ig,it),simTY1,indYT1,simTY2,indYT2)
      indYT(i2,i1,ig) = indYT1
      indYT(i2,nsim/2+i1,ig) = indYT2
      simY(i2,i1,ig) = simTY1*simPY2(i2,i1,ig)
      simY(i2,nsim/2+i1,ig) = simTY2*simPY2(i2,nsim/2+i1,ig)
   else
      indYT(i2,i1,ig) = 4
      indYT(i2,nsim/2+i1,ig) = 4 !same shock to guarantee that prob_unemp is matched
     simY(i2,i1,ig) = (1-unemp_dur)*simPY2(i2,i1,ig)+unemp_dur*min(yu1*simPY2(i2,i1,ig),yu1_cap)
     simY(i2,nsim/2+i1,ig) = (1-unemp_dur)*simPY2(i2,nsim/2+i1,ig)+unemp_dur*min(yu1*simPY2(i2,nsim/2+i1,ig),yu1_cap)
   end if    
   
end do !ig
end do !i2
end do  !i1 (nsim)

do ig=1,ng
    meanY(:,ig) = MATMUL(simY(:,:,ig),ones_nsim_1(:,1))/nsim
end do

!indices and coefficients for linear interpolation of income in the simulations
do i1=1,nsim
do t=1,tw+nper
do ig=1,ng
!the state variable "ginc" is only permanent income. We treat the transitory shock separately           
   simPY2(t,i1,ig) =  MAX(MIN(ginc(ninc,1),simPY2(t,i1,ig)),ginc(1,1))
   CALL ntoil(simPY2(t,i1,ig),ginc(:,1),ninc,i_inc)
   indy(t,i1,ig) = i_inc 
   tty(t,i1,ig) = (simPY2(t,i1,ig)-ginc(i_inc,1))/(ginc(i_inc+1,1)-ginc(i_inc,1))
   tty(t,i1,ig) = MAX(MIN(1.0,tty(t,i1,ig)),0.0)
end do
end do
end do

!POST-ROBOTIZATION
it=2

!cumulative densitiy function for the income shocks (with the top mass "empty" corresponding to the unemployment state)
do ig=1,ng
w1(1,1) = weig(1,1)*(1-p_unemp(ig,it))
do i7=2,n
   w1(i7,1) = w1(i7-1,1)+weig(i7,1)*(1-p_unemp(ig,it))
end do   
end do

do per=1,nper

do t=1,tw
do i1=1,nsim/2
   call random_number(rand_u) !before ig loop: same shocks for all groups
   call random_number(rand_u2) !before ig loop: same shocks for all groups
   DO ig=1,ng      
     if (per.eq.1) then
         if (indYT(t,i1,ig).eq.4) then
             simPY_1 = py_unemp
             simPY_2 = py_unemp
         else
            call simYP_f(rand_u2,weig(:,1),grid(:,1),n,smav(ig,it),simPY_1,simPY_2)
         end if
         simPY_R(per,t,i1,ig) = simPY_1+simPY(t,i1,ig)
         simPY_R(per,t,nsim/2+i1,ig) = simPY_2+simPY(t,nsim/2+i1,ig)
      else
         if (indYT_R(per-1,t,i1,ig).eq.4) then
             simPY_1 = py_unemp
             simPY_2 = py_unemp
         else
            call simYP_f(rand_u2,weig(:,1),grid(:,1),n,smav(ig,it),simPY_1,simPY_2)
         end if
         simPY_R(per,t,i1,ig) = simPY_1+simPY_R(per-1,t,i1,ig)
         simPY_R(per,t,nsim/2+i1,ig) = simPY_2+simPY_R(per-1,t,nsim/2+i1,ig)
      end if
      simPY2_R(per,t,i1,ig) = exp(simPY_R(per,t,i1,ig))*f_y(t+per,ig,it)
      simPY2_R(per,t,nsim/2+i1,ig) = exp(simPY_R(per,t,nsim/2+i1,ig))*f_y(t+per,ig,it)
    !for the transitory shocks we consider the unemployment and non-unemployment states separately
     if (rand_u<w1(3,1)) then
        call simY_f(rand_u,w1(:,1),grid(:,1),n+1,smay(ig,it),simTY1,indYT1,simTY2,indYT2)
        indYT_R(per,t,i1,ig) = indYT1
        indYT_R(per,t,nsim/2+i1,ig) = indYT2
        simY_R(per,t,i1,ig) = simTY1*simPY2_R(per,t,i1,ig)
        simY_R(per,t,nsim/2+i1,ig) = simTY2*simPY2_R(per,t,nsim/2+i1,ig)
     else
        indYT_R(per,t,i1,ig) = 4
        indYT_R(per,t,nsim/2+i1,ig) = 4 !same shock to guarantee that prob_unemp is matched
        simY_R(per,t,i1,ig) = (1-unemp_dur)*simPY2_R(per,t,i1,ig)+unemp_dur*min(yu1*simPY2_R(per,t,i1,ig),yu1_cap)
        simY_R(per,t,nsim/2+i1,ig) = (1-unemp_dur)*simPY2_R(per,t,nsim/2+i1,ig)+unemp_dur*min(yu1*simPY2_R(per,t,nsim/2+i1,ig),yu1_cap)
     end if 
   end do    
end do  !i1 (nsim)
end do !t

end do !per

!indices and coefficients for linear interpolation of income in the simulations
do per=1,nper
do t=1,tw
do i1=1,nsim
do ig=1,ng
   simPY2_R(per,t,i1,ig) =  MAX(MIN(ginc(ninc,1),simPY2_R(per,t,i1,ig)),ginc(1,1))
   CALL ntoil(simPY2_R(per,t,i1,ig),ginc(:,1),ninc,i_inc)
   indy_R(per,t,i1,ig) = i_inc 
   tty_R(per,t,i1,ig) = (simPY2_R(per,t,i1,ig)-ginc(i_inc,1))/(ginc(i_inc+1,1)-ginc(i_inc,1))
   tty_R(per,t,i1,ig) = MAX(MIN(1.0,tty_R(per,t,i1,ig)),0.0)
end do
end do
end do
end do

auxmyr(:,:,:) = simY_R(nper,:,:,:)
do ig=1,ng
do t=1,tw
    meanY_R(t,ig) = sum(auxmyr(t,:,ig))/nsim
end do
end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MAIN LOOPS: TYPES/INDUSTRIES and PRE/POST !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DO it=1,nt !PRE/POST
do ig=1,ng !INDUSTRY

write(*,*) ig, it

if (ig.eq.1) then
   rho = rho_L
   delta = delta_L
   psi = psi_L
else
   if (ig.eq.2) then
      rho = rho_M
      delta = delta_M
      psi = psi_M  
   else
      rho = rho_H
      delta = delta_H
      psi = psi_H
   end if
end if

theta = (1.0-rho)/(1.0-1.0/psi)
psi_1 = 1.0-1.0/psi
psi_2 = 1.0/psi_1 


!applicable when the CURRENT transitory income shock IS NOT an unemployment state
!note: i7 below is next-period's transitory income state, not the current one 
do i6=1,n
   do i8=1,n
      do i7=1,n
         nweig1(i6,i7,i8) = weig(i6,1)*weig(i7,1)*weig(i8,1)*(1-p_unemp(ig,it))*(1-p_lowr)
      end do
      nweig1(i6,n+1,i8) = weig(i6,1)*weig(i8,1)*p_unemp(ig,it)*(1-p_lowr)   
   end do
   do i7=1,n
      nweig1(i6,i7,n+1) = weig(i6,1)*weig(i7,1)*p_lowr*(1-p_unemp(ig,it))
   end do
   nweig1(i6,n+1,n+1) = weig(i6,1)*p_unemp(ig,it)*p_lowr     
end do
!applicable when the CURRENT transitory income shock IS an unemployment state
do i8=1,n
   do i7=1,n
      nweig1b(i7,i8) = weig(i7,1)*weig(i8,1)*(1-p_unemp(ig,it))*(1-p_lowr)
   end do
   nweig1b(n+1,i8) = weig(i8,1)*p_unemp(ig,it)*(1-p_lowr)   
end do
do i7=1,n
   nweig1b(i7,n+1) = weig(i7,1)*p_lowr*(1-p_unemp(ig,it))
end do
nweig1b(n+1,n+1) = p_unemp(ig,it)*p_lowr     


!!!!!!!!!!!!!!!!!!!!!!!
! LABOR INCOME SHOCKS !
!!!!!!!!!!!!!!!!!!!!!!!
!unemployment shock is captured separately
do i1=1,n
   yh(i1,1) = EXP(grid(i1,1)*smay(ig,it))
end do

do i1=1,n
   grid2(:,1) = grid(i1,1)*corr_v+grid(:,1)*ones_n_1(:,1)*(1-corr_v**2)**(0.5)
   yp(:,i1) = exp(grid2(:,1)*smav(ig,it))
end do
!for simplicity, set income shock equal to mean value in return disaster state (reduces correlation/risk sightly)
yp(:,n+1) = EXP(grid((n+1)/2,1)*smav(ig,it)) 

!!!!!!!!!!!!!!!!!!!!!!!
! LABOR INCOME VALUES !
!!!!!!!!!!!!!!!!!!!!!!!
!when current state IS NOT an unemployment state
do i2=1,tr-tb-1
do iy=1,ninc
do i6=1,n
    do i8=1,n+1 !return 
       inc_t1(iy,i6,i8,i2) = yp(i6,i8)*ginc(iy,1)*f_y(i2+1,ig,it)/f_y(i2,ig,it)
    end do
end do !i6
end do !iy
end do !i2
do iy=1,ninc
   inc_t1(iy,:,:,tr-tb) = ret_fac0(ig,it)*ginc(iy,1)
end do !iy

do i2=1,tr-tb
do iy=1,ninc
do i6=1,n
   do i8=1,n+1 !return 
      inc_t1(iy,i6,i8,i2) =  MAX(MIN(ginc(ninc,1),inc_t1(iy,i6,i8,i2)),ginc(1,1))         
      CALL ntoil(inc_t1(iy,i6,i8,i2),ginc(:,1),ninc,i_inc)
      inc_it1(iy,i6,i8,i2) = i_inc 
      tty_t1(iy,i6,i8,i2) = (inc_t1(iy,i6,i8,i2)-ginc(i_inc,1))/(ginc(i_inc+1,1)-ginc(i_inc,1))
      tty_t1(iy,i6,i8,i2) = MAX(MIN(1.0,tty_t1(iy,i6,i8,i2)),0.0)
  end do
end do !i6
end do !iy
end do !i2

!when current state IS an unemployment state
do i2=1,tr-tb-1
do iy=1,ninc
   inc_t1_u(iy,i2) = exp(py_unemp)*ginc(iy,1)*f_y(i2+1,ig,it)/f_y(i2,ig,it)
end do !iy
end do !i2
do iy=1,ninc
   inc_t1_u(iy,tr-tb) = ret_fac0(ig,it)*ginc(iy,1)
end do !iy


do i2=1,tr-tb
do iy=1,ninc
   inc_t1_u(iy,i2) =  MAX(MIN(ginc(ninc,1),inc_t1_u(iy,i2)),ginc(1,1))         
   CALL ntoil(inc_t1_u(iy,i2),ginc(:,1),ninc,i_inc)
   inc_it1_u(iy,i2) = i_inc 
   tty_t1_u(iy,i2) = (inc_t1_u(iy,i2)-ginc(i_inc,1))/(ginc(i_inc+1,1)-ginc(i_inc,1))
   tty_t1_u(iy,i2) = MAX(MIN(1.0,tty_t1_u(iy,i2)),0.0)
end do !iy
end do !i2

!!!!!!!!!!!!!!!!!!!
! TERMINAL PERIOD !
!!!!!!!!!!!!!!!!!!!
do i1=1,nlw
do iy=1,ninc
   C(i1,tn,iy,:,ig,it) = glw(i1,1)+ginc(iy,1)  !same for all values of the transitory shock
end do
end do
do i1=1,nlw
do iy=1,ninc
do i7=1,n+1    
   V(i1,tn,iy,i7,ig,it) = C(i1,tn,iy,i7,ig,it)*((1.0-delta)**(psi/(psi-1.0)))
end do
end do
end do
A(:,tn,:,:,ig,it) = 0.0


do iy=1,ninc
   call spline(glw(:,1),V(:,tn,iy,1,ig,it),nlw,1.0,secd(:,iy,1,ig,it))
   do i7=2,n+1    
      secd(:,iy,i7,ig,it) = secd(:,iy,1,ig,it)
   end do
end do


!!!!!!!!!!!!!!!!!!!!!!
! RETIREMENT PERIODS !
!!!!!!!!!!!!!!!!!!!!!!
DO i1= 1,td-tr
   t= tn-i1
   write(*,*) t, ig, it
   do iy=1,ninc
   survprob = min(max(survprob0(t,1)+survprob1(t,1)*ginc(iy,1),0.00001),0.99999)     
   do i3=1,nlw
      call cgrid(i3,C(1,t+1,iy,1,ig,it),nlw,C(:,t,iy,1,ig,it),glw,nc,gc)
      do i4=1,nc
         caux = gc(i4,1)
         u=(1.0-delta)*(caux**psi_1) !1-1/psi
         sav = glw(i3,1)+(1-Fy)*ginc(iy,1)-gc(i4,1)-F
         do i5=1,na
            if (i5.eq.na) then
               sav = sav+F+Fy*ginc(iy,1)
            end if
call EVret(n+1,riskret(i5,:),sav,glw,V(:,t+1,iy,1,ig,it),secd(:,iy,1,ig,it),nlw,nweig2,survprob,rho,auxVV)
            auxV(i5,i4) = (u+delta*(auxVV**(1.0/theta)))**psi_2    !1/(1-1/psi)
         end do    !i5
      end do   !i4
      vec_V = RESHAPE((auxV),(/na*nc,1/))
      V(i3,t,iy,1,ig,it) = MAXVAL(vec_V(:,1))
      pt = MAXLOC(vec_V(:,1))
      aux2 = FLOOR((REAL(pt(1))-1)/REAL(na))
      C(i3,t,iy,1,ig,it) = gc(aux2+1,1)
      A(i3,t,iy,1,ig,it) = ga(pt(1)-na*aux2,1)
    end do    !i3

   call spline(glw(:,1),V(:,t,iy,1,ig,it),nlw,1.0,secd(:,iy,1,ig,it))
   end do !iy
   
   do i7=2,n+1    
      V(:,t,:,i7,ig,it) = V(:,t,:,1,ig,it)
      C(:,t,:,i7,ig,it) = C(:,t,:,1,ig,it)
      A(:,t,:,i7,ig,it) = A(:,t,:,1,ig,it)
      secd(:,:,i7,ig,it) = secd(:,:,1,ig,it)
   end do

   
end do    !i1


!!!!!!!!!!!!!!!!!
! OTHER PERIODS !
!!!!!!!!!!!!!!!!!
DO i1= 1,tr-tb
   t= tr-tb-i1+1
   write(*,*) t, ig, it
   do iy=1,ninc !permanent income
   survprob = min(max(survprob0(t,1)+survprob1(t,1)*ginc(iy,1),0.00001),0.99999)  
   do i7=1,n+1   !transitory income shock 
      if (i7>n) then
         income = (1-unemp_dur)*ginc(iy,1)+unemp_dur*min(yu1*ginc(iy,1),yu1_cap)
      else
         income =  ginc(iy,1)*yh(i7,1)           
      end if 
   do i3=1,nlw
      call cgrid(i3,C(1,t+1,iy,i7,ig,it),nlw,C(:,t,iy,i7,ig,it),glw,nc,gc)
      do i4=1,nc
         caux = gc(i4,1)
         u=(1.0-delta)*(caux**psi_1)       !1-1/psi
         sav = glw(i3,1)+(1-Fy)*income-gc(i4,1)-F   
          do i5=1,na
            if (i5.eq.na) then
               sav=sav+F+Fy*ginc(iy,1)
            else
              if (sav<0) then
                 gc(i4,1)= gc(i4,1)+sav !adding a negative number  
                 sav = 0.0
              end if                 
            end if 
if (i7>n) then
    call EVu(n+1,n+1,riskret(i5,:),sav,inc_it1_u(iy,t),tty_t1_u(iy,t),glw,nlw,ninc, &
    V(:,t+1,:,:,ig,it),secd(:,:,:,ig,it),nweig1b,survprob,rho,auxVV)    
else    
    call EV(n,n+1,n+1,riskret(i5,:),sav,inc_it1(iy,:,:,t),tty_t1(iy,:,:,t),glw,nlw,ninc, &
    V(:,t+1,:,:,ig,it),secd(:,:,:,ig,it),nweig1,survprob,rho,auxVV)
end if
             auxV(i5,i4) = (u+delta*(auxVV**(1.0/theta)))**psi_2    !1/(1-1/psi)     
         end do    !i5
      end do   !i4     
      vec_V = RESHAPE((auxV),(/na*nc,1/))
      V(i3,t,iy,i7,ig,it) = MAXVAL(vec_V(:,1))
      pt = MAXLOC(vec_V(:,1))
      aux2 = FLOOR((REAL(pt(1))-1)/REAL(na))
      C(i3,t,iy,i7,ig,it) = gc(aux2+1,1)
      A(i3,t,iy,i7,ig,it) = ga(pt(1)-na*aux2,1)   
      
   end do    !i3
 
!   call spline(glw(:,1),V(:,t,iy,1,ig,it),nlw,1.0,secd(:,iy,1,ig,it))
   call spline(glw(:,1),V(:,t,iy,i7,ig,it),nlw,1.0,secd(:,iy,i7,ig,it))
   end do !i7
 end do !iy
    
end do    !i1

end do !ig
end do !it


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATIONS !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do ig=1,ng !INDUSTRY

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATIONS: PRE-ROBOTIZATION !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
it=1 

   write(*,*) ig, it

simW(:,:)=0.0
simP(:,:)=0.0
simW(1,:) = W0
do t=1,tw

do i1=1,nsim
   simW_Y(t,i1) = simW(t,i1)/simY(t,i1,ig)
   cash = max(simW(t,i1),glw(1,1))
   income = simY(t,i1,ig)
   
   call sim_dec(nlw,ninc,glw,lglw,C(:,t,:,indYT(t,i1,ig),ig,it),A(:,t,:,indYT(t,i1,ig),ig,it),cash,income,tty(t,i1,ig),indy(t,i1,ig),F,Fy, &
       simC(t,i1),simA(t,i1),simS(t,i1),simB(t,i1),simBpart(t,i1),simP(t,i1))
   if (t<tw) then
      simW(t+1,i1) = simB(t,i1)*r+simS(t,i1)*simR(t,i1)
   end if

 
end do  !i1 (nsim)

end do  !t 

meanC(:,ig) = MATMUL(simC(:,:),ones_nsim_1(:,1))/nsim
meanW(:,ig) = MATMUL(simW(:,:),ones_nsim_1(:,1))/nsim
meanS(:,ig) = MATMUL(simS(:,:),ones_nsim_1(:,1))/nsim
meanB(:,ig) = MATMUL(simB(:,:),ones_nsim_1(:,1))/nsim
meanBpart(:,ig) = MATMUL(simBpart(:,:),ones_nsim_1(:,1))/nsim
meanP(:,ig) = MATMUL(simP(:,:),ones_nsim_1(:,1))/nsim
meanWY(:,ig) = MATMUL(simW_Y(:,:),ones_nsim_1(:,1))/nsim
!meanalpha(:,1) = MATMUL(simA(:,:),ones_nsim_1(:,1))/nsim

do t=1,tw
   meanalpha(t,ig) = 0.0
   partcount=0.0
   do i1=1,nsim
      if (simP(t,i1)>0) then
          meanalpha(t,ig) = meanalpha(t,ig)+simA(t,i1)
          partcount = partcount+1.0
       end if    
   end do  !i1 (nsim)
   if (partcount>0.0) then
       meanalpha(t,ig) = meanalpha(t,ig)/partcount
   else    
       meanalpha(t,ig) = 0.0 !irrelevant (just to avoid the NaN above)
   end if    
end do  !t 

mWY(ig,1) = 0.0
mW(ig,1) = 0.0
mP(ig,1) = 0.0
malpha(ig,1) = 0.0
do t=1,tw
   mW(ig,1) = meanW(t,ig)*dweig(t,1)/100 + mW(ig,1)
   mWY(ig,1) = meanWY(t,ig)*dweig(t,1)/100 + mWY(ig,1)
   mP(ig,1) = meanP(t,ig)*dweig(t,1)/100 + mP(ig,1)
   malpha(ig,1) = meanalpha(t,ig)*dweig(t,1)/100 + malpha(ig,1)
end do

!ng+1 slot is not re-set, it "adds"/averages with each ig loop (hence divided by 300 instead of 100)
do t=1,tw
   mW(ng+1,1) = gs(ig,1)*meanW(t,ig)*dweig(t,1)/100 + mW(ng+1,1)
   mWY(ng+1,1) = gs(ig,1)*meanWY(t,ig)*dweig(t,1)/100 + mWY(ng+1,1)
   mP(ng+1,1) = gs(ig,1)*meanP(t,ig)*dweig(t,1)/100 + mP(ng+1,1)
   malpha(ng+1,1) = gs(ig,1)*meanalpha(t,ig)*dweig(t,1)/100 + malpha(ng+1,1)
end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATIONS: ROBOTIZATION SHOCK !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

it=2
   write(*,*) ig, it
do t=1,tw  !DATE AT WHICH THE SHOCK HITS

do per=1,nper !number of periods after the shock
   t2=t+per 

do i1=1,nsim
   if (per.eq.1) then
      cash = max(simW(t,i1),glw(1,1))
      simW_Y_aux = simW(t,i1)/simY_R(per,t,i1,ig)
   else
      cash = max(simW_aux(per,i1),glw(1,1))
      simW_Y_aux = simW_aux(per,i1)/simY_R(per,t,i1,ig)
   end if   
   income = simY_R(per,t,i1,ig)
   call sim_dec(nlw,ninc,glw,lglw,C(:,t,:,indYT_R(per,t,i1,ig),ig,it),A(:,t,:,indYT_R(per,t,i1,ig),ig,it),cash,income,tty_R(per,t,i1,ig),indy_R(per,t,i1,ig),F,Fy, &
       simC_aux,simA_aux,simS_aux,simB_aux,simBpart_aux,simP_aux)
   if (per<nper) then
      simW_aux(per+1,i1) = simB_aux*r+simS_aux*(simR(t+per,i1)+Radj)
   end if

 
end do  !i1 (nsim)

end do ! per

simW_R(t,:) = simW_aux(nper,:)

end do  !t 


meanW_R(:,ig) = MATMUL(simW_R(:,:),ones_nsim_1(:,1))/nsim

!Note: we are comparing with the same person in age t, so the demograpic weights have to be the same (same individuals after all)
mW_R(ig,1) = 0.0
do t=1,tw
   mW_R(ig,1) = meanW_R(t,ig)*dweig(t,1)/100 + mW_R(ig,1)
end do

do t=1,tw
   mW_R(ng+1,1) = gs(ig,1)*meanW_R(t,ig)*dweig(t,1)/100 + mW_R(ng+1,1)
end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATIONS: ROBOTIZATION SHOCK WITH OLD ALPHA POLICY FUNCTIONS !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

it=2
   write(*,*) ig, it
do t=1,tw  !DATE AT WHICH THE SHOCK HITS

do per=1,nper !number of periods after the shock
   t2=t+per 

do i1=1,nsim
   if (per.eq.1) then
      cash = max(simW(t,i1),glw(1,1))
      simW_Y_aux = simW(t,i1)/simY_R(per,t,i1,ig)
   else
      cash = max(simW_aux(per,i1),glw(1,1))
      simW_Y_aux = simW_aux(per,i1)/simY_R(per,t,i1,ig)
   end if   
   income = simY_R(per,t,i1,ig)
   call sim_dec(nlw,ninc,glw,lglw,C(:,t,:,indYT_R(per,t,i1,ig),ig,it),A(:,t,:,indYT_R(per,t,i1,ig),ig,1),cash,income,tty_R(per,t,i1,ig),indy_R(per,t,i1,ig),F,Fy, &
       simC_aux,simA_aux,simS_aux,simB_aux,simBpart_aux,simP_aux)
   if (per<nper) then
      simW_aux(per+1,i1) = simB_aux*r+simS_aux*(simR(t+per,i1)+Radj)
   end if

end do  !i1 (nsim)

end do ! per

simW_R(t,:) = simW_aux(nper,:)

end do  !t 


meanW_R2(:,ig) = MATMUL(simW_R(:,:),ones_nsim_1(:,1))/nsim

mW_R2(ig,1) = 0.0
do t=1,tw
   mW_R2(ig,1) = meanW_R2(t,ig)*dweig(t,1)/100 + mW_R2(ig,1)
end do

do t=1,tw
   mW_R2(ng+1,1) = gs(ig,1)*meanW_R2(t,ig)*dweig(t,1)/100 + mW_R2(ng+1,1)
end do

end do !ig

OPEN(UNIT=25,FILE='Sum',STATUS='replace',ACTION='write')
WRITE(25,*) mWY(:,1) 
WRITE(25,*) mP(:,1)
WRITE(25,*) malpha(:,1)
WRITE(25,*) mW(:,1) 
WRITE(25,*) mW_R(:,1)
WRITE(25,*) mW_R2(:,1)


CLOSE(UNIT=25)


end program

    